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Abstract 

Pituitary adenomas (PAs) are neoplasms that may cause a variety of neurological and endocrine effects. Although known 
causal contributors include heredity, hormonal influence and somatic mutations, the pathophysiologic mechanisms driving 
tumorigenesis and invasion of sporadic PAs remain unknown. We hypothesized that alterations in DNA methylation are 
associated with PA invasion and histopathology subtype, and that genome-scale methylation analysis may complement 
current classification methods for sporadic PAs. Twenty-four surgically-resected sporadic PAs with varying histopathological 
subtypes were assigned dichotomized Knosp invasion scores and examined using genome-wide DNA methylation profiling 
and RNA sequencing. PA samples clustered into subgroups according to functional status. Compared with hormonally- 
active PAs, nonfunctional PAs exhibited global DNA hypermethylation (mean beta-value 0.47 versus 0.42, P = 0.005); the 
most significant site of differential DNA methylation was within the promoter region of the potassium voltage-gated 
channel KCNAB2 (FDR = 5.11 x10~ 10 ). Pathway analysis of promoter-associated CpGs showed that nonfunctional PAs are 
potentially associated with the ion-channel activity signal pathway. DNA hypermethylation tended to be negatively 
correlated with gene expression. DNA methylation analysis may be used to identify candidate genes involved in PA function 
and may potentially complement current standard immunostaining classification in sporadic PAs. DNA hypermethylation of 
KCNAB2 and downstream ion-channel activity signal pathways may contribute to the endocrine-inactive status of 
nonfunctional PAs. 
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Introduction 

Pituitary adenomas (PAs) are typically benign monoclonal 
neoplasms with an overall prevalence of 16.7% (14.4% in autopsy 
studies and 22.5% in imaging studies) in the general population. 
The majority of PAs, however, are small and nonfunctional 
tumors, and only 0.16-0.20% of them are macroadenomas S 
10 mm in diameter [1,2]. Although histologically benign, many 
PAs may cause significant morbidity due to their anatomical 
location, often resulting in tumor mass effect and neurological 



symptoms in addition to causing hormonal over-secretion or 
hypopituitarism. Invasion into surrounding anatomical structures 
(e.g. cavernous sinus invasion) remains a major barrier to 
achieving long-term tumor and disease control, especially in cases 
of functional PAs resulting in malignant endocrinopathies such as 
Cushing's disease or acromegaly [3]. With regard to the 
phenotype of invasion, PAs are often classified according to the 
Knosp grading system, in which one of five grades is assigned 
based on the degree of cavernous sinus invasion and relationship 
to the internal carotid artery [4]. 
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Previous studies have shown that genetic mutations play an 
important role in the tumorigenesis of familial PAs, particularly 
those with heritable mutations in the multiple endocrine neoplasia 
I (MEN1) and aryl hydrocarbon receptor interacting protein (AIP) 
genes [5]. In sporadic PAs, however, it has been suggested that 
alterations in epigenetic regulation, and particularly DNA 
methylation, may play a particularly prominent role in PA 
tumorigenesis and invasion, likely via loss or reduced expression of 
tumor suppressor genes (TSGs) [6] . Many known TSGs have been 
shown to harbor C-phosphate-G (CpG) island hypermethylation 
in PAs, which is associated with and frequently results in silencing 
of TSGs. For example, DNA methylation within promoter and 
exon 1 regions of the cyclin-dependent kinase inhibitor {pi 6/ 
CDKN2A) gene was found at a high frequency in NFAs [7,8], and 
DNA methylation of the CpG island like cell-cycle regulatory 
genes of growth arrest and DNA-damage-inducible, gamma 
[GADD45G) [9] and apoptosis gene of rhomboid domain 
containing 3 (RHBDD3) have been associated with PA evolution 
[10]. Other studies have shown that tumor specific epigenetic 
silencing of cadherin 13, H-cadherin (CDH13) and cadherin 1, 
typel, E-cadherin (CDH1), alone or in combination, were involved 
in PA development and invasion [11], and CpG hypermethyla- 
tion-mediated glutathione S-transferase pi gene (GSTP1) inactiva- 
tion was a common finding in PAs potentially contributing to their 
invasive behavior [12]. Few prior studies, however, have utilized 
genome-scale approaches according to functional PA subtypes and 
phenotypical invasion status to identify candidate genes involved 
in these processes [13,14]. 

Although several targeted genes have been associated with PA 
tumorigenesis and progression, primary whole genome-scale 
epigenetic alterations remain largely unknown, especially as they 
pertain to invasion and histopathological PA subtype classification. 
In the current study, we utilized genome-scale DNA methylation 
technology to identify DNA methylation alterations between 
invasive and noninvasive PAs subtypes, and across varying 
functional classes of PAs. 

Materials and Methods 

Patients and tissue specimens 

This study was conducted according to the Helsinki human 
subject doctrine and was approved by the Institutional Review 
Boards of the Keck School of Medicine of the University of 
Southern California (Los Angeles, USA). Written informed 
consent was signed and obtained from 23 participants and one 
guardian on behalf of a minor enrolled for tissue specimen 
collection and subsequent analysis. Twenty-four patients with 
surgically-resected PAs from the University of Southern California 
(USC) Keck Hospital and Los Angeles County + USC Medical 
Center were included in this retrospective study. 

Of the 24 PAs, six (25%) were functional and 18 (75%) were 
nonfunctional. More specifically, these included 17 (70.8%) 
nonfunctional adenomas, 5 (20.8%) somatotroph adenomas, 1 
(4.2%) corticotroph adenoma, and 1 (4.2%) silent corticotroph 
adenoma. All PAs analyzed were macroadenomas > 10 mm in 
diameter, and were diagnosed based on laboratory evaluation and 
Magnetic Resonance Imaging (MRI), with an average tumor 
diameter of 23.55 TRx 19.36 CC mm (Table 1). Knosp invasion 
scoring on MRI was performed by a staff neuro-radiologist who 
was blinded to the genomics analysis. For the purposes of this 
study, PAs with Knosp scores of 0-1 were classified as noninvasive, 
and those with Knosp scores of 2-4 were classified as invasive. 
Equal numbers of invasive and noninvasive PAs (12 each) were 
selected. 



Following transsphenoidal tumor resection, PA specimens were 
fresh-frozen in liquid nitrogen and stored at — 80°C until used for 
DNA and RNA extraction. Genomic DNA purification was 
processed using DNeasy Blood & Tissue kit (Qiagen) according to 
the standard protocols of the manufacturer. Total RNA was 
purified using RNeasy Plus Universal Kits (Qiagen). Tissue 
samples were first cut into small pieces with no more than 
50 mg on dry ice then moved promptly into 1.5 ml precooled 
tubes. Liquid nitrogen was then added into the tube, and tissue 
was grinded with a grinding rod before adding 0.9 ml of QIAzol 
Lysis Reagent. The tubes were vortexed for 60 seconds and left at 
room temperature for 5 minutes before following the standard 
RNeasy Plus protocol for RNA purification. 

Genome-scale DNA methylation array 

One microgram of bisulphite converted DNA (using the Zymo 
EZ DNA Methylation kit, Zymo Research, Irvine, CA) per sample 
was used for genome-scale DNA methylation profiling, and was 
performed in the USC Epigenome Center using the Illumina 
Infinium HumanMethylation450 (HM450) Beadchip platform, 
which interrogates 482,421 DNA methylation sites and covers 
99% of Reference Sequence of National Center for Biotechnology 
Information (NCBI-RefSeq) genes, with an average of 17 CpG 
sites per gene region distributed across the promoter, TSS1500 
(1500 bp within transcription start site), TSS200 (200 bp within 
transcription start site), 5'UTR, IstExon (the first exon), gene 
body, 3'UTR. It covers 96% of CpG islands, with additional 
coverage in island shores and the regions flanking them. The DNA 
methylation data was submitted to Gene Expression Omnibus 
(GEO) with access number of GSE54415. Beta values for all 
HM450 probes were calculated after background correction and 
normalization as previously described [15]. Probes (11,648) 
located on chromosomes X and Y were excluded. Probes targeting 
non-CpG sites (3,091) and probes targeting 65 known single 
nucleotide polymorphisms (SNPs) were also excluded. An addi- 
tional set of 86,560 probes with common SNPs at the CpG site 
with a minor allele frequency (MAF) greater than 1% (identified 
using NCBI dbSNP builds 128), or with common SNPs within 
10 bp from CpG site, or within 15 bp from the CpG site lying 
entirely within a repeat region (from RepeatMaster and Tandem 
Repeat Finder databases) were excluded. Finally, 591 probes with 
detection P-values >0.01 were removed. Following these filtering 
steps, a total of 383,718 probes were utilized for analyses (Figure A 
in File SI). We determined the distribution of each probe by CpG 
density (island, shore, shelf, etc) and by location (promoter, gene 
body, UTR, etc.) (Figure A and B in File S2). A secondary analysis 
of probes with beta values <0.1 across 77.8% (14/18 samples) of 
nonfunctional adenomas were found to be less variable and 
therefore could conceal the potential variance. These probes were 
excluded in order to investigate the underlying difference between 
invasive and non-invasive nonfunctional adenomas. Remaining 
subsets of 884 probes were selected for the in-depth analysis 
(Figure B in File SI). Finally, probes with Abeta >± 0.3 were used 
in differentiating various PA classes based on phenotype. 

Hierarchical Clustering (HCL) was performed using the 
Microarray Software Suite of MultiExperiment Viewer (MeV) 
V4.8 [16] using probes with the top 1%, 2% and 5% of standard 
deviation (SD) across the set of 383,718 "global" probes in 24 PA 
samples respectively. Two-way clustering was performed using 
Pearson's correlation and average linkage for both the gene tree 
and sample tree (Figure 1 with top 2% SD). Differential DNA 
methylation analyses were performed between: 1) Invasive versus 
noninvasive PAs (12 samples each), 2) FA subjects (n = 6) versus 
NFA subjects (n= 18), 3) Invasive NFAs (n= 10) versus noninva- 
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sive NFAs (n = 8), and 4) PAs and NFAs stratified according to the 
five different Knosp invasion grades (0-4). Two tailed t-test 
analyses were performed between FAs versus NFAs and noninva- 
sion versus invasion PA. Kruskal-Wallis test was subsequently used 
to analyze DNA methylation differences among 18 NFA subjects 
with five different Knosp grades (Table 1). All of the analysis were 
executed in globe probes with the proportion of false significant 
genes did not exceed 0.01. 

MethyLight DNA methylation analysis 

Validation of findings from DNA methylation profiling in 
KCNAB2, which spans the CpG site interrogated by the Illumina 
HM450 array at probe cgl8 192083, was performed using 
MethyLight real-time PCR assays as previously described [17]. 
One microgram of genomic DNA from each of the 24 PA samples 
and the CpG methyltransferase (M.&T) treated reference DNA 
sample was converted with bisulfite using the Zymo EZ DNA 
Methylation kit (Zymo Research, Irvine, CA), as specified by the 
manufacturer. A control reaction targeting ALU repetitive 
elements [18] was used to normalize the input bisulfite-DNA 



amounts. The real-time PCR reactions were obtained from Life 
Technologies, Carlsbad, CA. The KC.NAB2 MethyLight primer 
and probe sequences are: forward: 5' GGT TTT TTT ATT TGG 
GTT ACG CG-3'; reverse: 5'-CGC AAA ACT AAA AAA CCT 
AAC GC-3'; probe: J'-6FAM-CCG TAA AAT ATC GAA ACG 
TAA CC-MGBNFQ;3', in which 6FAM is the fluorophore and 
MGBNFQ, refers to a minor-groove binding non-fluorescent 
quencher. The ALU control reaction primer and probe sequences 
were previously described [18]. The data were reported as a 
Percent of Methylated Reference (PMR), calculated as: PMR = 
{{{KCNAB2/ALU) S!SXap ^/{KCNAB2/ALU] M . Sss i)) x 100. 

Using HM450 manifest (downloaded from gene expression 
omnibus, GPL13534) 86,772 promoter-associated CpG probes 
were selected for the pathway analysis using gene set enrichment 
analysis (GSEA) software [19]. A total of 1,454 gene sets of 
Human Gene Ontology were exported from the molecular 
signatures database (MSigDB) version 3.1, and after filtering, 
469 of the gene sets were either larger than 500 or smaller than 25 
gene set size allowed. Additional statistical analysis was performed 
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using R software (http://www.r-project.org/) and GraphPad 
Prism Software (GraphPad Software Inc.). 

Whole-genome gene expression sequencing 

One microgram of total RNA from each of the 24 samples was 
assessed by Agilent Bioanalyzer to examine quality, and subse- 
quently used for library preparation. Library quantization was 
performed by Kapa Biosystems qPCR assay. RNA sequencing was 
performed in the USC Epigenome Center using Illumina Hi-seq 
2000 platform to produce 100 bp paired-end reads, by indexing 
four samples per lane. The Fastq files could be downloaded from 
Sequence Read Archive (SRA) with accession number 
SRP035646. Raw sequence data were first examined by FastQC 
as a quality control measure. We next used Tophat and Cuffdiff to 
examine genes with differential expression, based on the fold 
change of fragments of per kilo base of transcript sequence per 
millions base pairs (FPKM). 

Results 

Genome-wide DNA methylation in invasive and 
noninvasive PAs 

In this pilot study, twenty-four surgically-resected sporadic PAs 
with varying histopathological subtypes and Knosp invasion scores 
were characterized to determine if genome-wide and gene-specific 
methylation levels can serve as biomarkers of PA subtype and 
invasion grading. We calculated the mean DNA methylation beta 
value for the entire filtered probe set (383,7 1 8 probes) as a measure 
for global DNA methylation in noninvasive and invasive PAs 
(Figure 2A), across different tumor grades by Knosp categorization 
(Figure 2B), and in FAs and NFAs (Figure 2C). Based on these 
calculations, no significant differences in global DNA methylation 
levels were observed between the invasive and noninvasive groups 
(Figure 2A). Furthermore, Kruskal-Wallis test analysis among PAs 
with five different Knosp grades using the global filtered probe set 
also did not show significant DNA methylation differences as a 
function of tumor grade (Figure 2B), even though we relaxed the 
threshold to FDR <0.05. Within the set of 18 NFA subjects, 
similar (nonsignificant) findings between invasive and noninvasive 
PAs were also seen (Figure 2D). These data suggest that global 
DNA methylation profiles may not independently predict clinically 
significant differences in the invasive phenotype of PAs. 

We further refined the probe set for unsupervised analyses to 
remove non-variably methylated probes. After filtering, 884 
probes were selected, of which 355 were enhancer-associated, 20 
were promoter-associated, and 4 were both promoter- and 
enhancer-associated based on the Illumina HM450 manifest for 
probes (Figure B in File SI). We identified 34 CpGs (localized to 
1 7 genes) that were independendy associated with enhancers and 
were hypomethylated in invasive NFAs compared to noninvasive 
NFAs (FDR <0.01) (Table SI). These differentially-methylated 
genes included the fms-related tyrosine kinase 1 (FLT1) and slit 
homolog 3 (SLIT3), which are two important genes known to be 
associated with cell motility and invasion. Although we attempted 
to validate methylation status in these two genes using Methy- 
Light, we were not able to because the probes lie in CpG-depleted 
regions of the genes. No DNA methylation differences in promoter 
regions were noted between invasive and noninvasive NFAs. 

Nonfunctional PAs are globally hypermethylated 
compared to somatotroph PAs 

The group of 24 PA samples independently clustered into two 
subgroups of nonfunctional and functional (all somatotroph) PAs 
following HCL with top 2% of SD (Figure 1). All somatotroph 



adenomas were exclusively clustered from NFAs (and one 
additional functional corticotroph adenoma) (Figure 1), and 
further clustering of the top 1% (File S3) and 5% (File S4) of SD 
validated this finding. The mean global beta values from FA (six 
subjects) and NFA (18 subjects) tumors showed that NFAs were 
significandy DNA hypermethylated compared to FAs (Figure 2C). 
However, NFAs and FAs showed similar DNA methylation trends 
from gene promoters through gene bodies and UTR regions 
(Figure 2E and F). These findings suggest that DNA methylation 
analysis may provide insight into hormonal activity in functional 
versus nonfunctional PAs, and possibly complement current PA 
histopathological classification systems. 

In order to interrogate which subsets of genes were involved in 
distinguishing specialization of nonfunctional PAs, we filtered the 
global panel of 383,718 HM450 probes to include those targeted 
in RefSeq genes for a more thorough investigation. Using a two- 
tailed t-test, we identified 3,027 CpG sites with significant DNA 
methylation differences between NFAs and FAs (FDR <0.01). In 
order to display the methylation status of the significant CpG sites, 
we graphed DNA methylation levels of these probes stratified by 
location relative to gene region (Figure 3A) and as a function of 
CpG density (Figure 3B). We further restricted the set of 3,072 
probes to 360 CpG sites associated with promoter regions, and 
found substantial DNA hypomethylation in FA cases compared to 
NFA cases (simple t-test, P<0.01) (Figure 3C). Furthermore, CpG 
sites were stratified by density (Figure 3D) as well. The most 
significant gene with regard to differential DNA methylation was 
KCNAB2, a member of the potassium voltage-gated channel, 
shaker-related subfamily (FDR = 5. 1 1 x 10" lu ), which was hyper- 
methylated in the NFA group. This result was subsequendy 
validated using MethyLight real-time PCR assays (Figure 3E). 
Furthermore, we found that among the top 30 most significandy 
differentiated genes, 4 {NME9, PSEN2, MDN1 and PCSK6) were 
hypomethylated in NFAs compared to FAs (Table 2). Interesting- 
ly, we also corroborated the recent finding discovered using 
Illumina Infinium Methylation 27K Arrays that the transcription 
factor AP-2 epsilon (TFAP2E) and 13 other genes (Table S2) were 
hypermethylated in NFAs [13] (Table 2). Epigenetic alteration of 
the TFAP2E gene has also been commonly reported in selected 
human cancers, including prostate cancer and colorectal carcino- 
ma [20,21]. In addition, Wnt signaling plays an important role in 
adenoma development and tumorigenesis, including PAs [22]. 
Whether hypermethylation of TFAP2E could down-regulate Wnt 
signaling, however, was still unknown. In order to further 
understand the detailed DNA methylation status of KCNAB2, we 
performed hierarchical clustering of beta values from KGNAB2- 
targeted probes on the HM450 array using the top 50% highest 
standard deviation. In nonfunctional PAs, KCNAB2 was generally 
hypermethylated across the promoter region (TSS1500 and 
TSS200), 5 'UTR, the first exon, and gene body. This was 
substantially different from the coverage profile observed in 
functional PAs (File S5). Since KCMAB2 is expressed abundandy in 
the nervous system and T lymphocytes, and appears to play an 
important role in K + channel activation [23] and subsequent 
hormone secretion [24], we further extended our investigation of 
this pathway. 

Methylation of ion channel activity signal pathway genes 
may be associated with PA functional status 

Because gene promoter DNA hypermethylation is known to 
play an important role in tumorigenesis [25] and contributes to 
loss of gene function and silencing [26], we selected promoter- 
associated probes according to the Illumina HM450 manifest for 
further pathway analyses. We found that four gene sets were 
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Figure 2. Global CpG DNA methylation profiles in PAs. A, Comparison of mean DNA methylation levels across 383,718 HM450 probes in 
noninvasive (Knosp grade 0, 1) and invasive (Knosp grade 2, 3, 4) PAs. B, Mean DNA methylation levels in PA cases across the five different Knosp 
grades (0-4) among the set of 383,718 HM450 probes. C, Comparison of mean DNA methylation levels across 383,718 HM450 probes in FAs and 
NFAs. P-value was calculated via simple t-test. D, Comparison of mean DNA methylation levels across 383,718 HIV1450 probes in only NFAs. P-value 
was calculated via simple t-test. All of the above DNA methylation beta values are shown as mean ± SEM in different groups. E and F, DNA 
methylation levels in gene regions (TSS1500, TSS200, 5'UTR, IstExon, gene body and 3'UTR) for NFAs (panel E) and FAs (Panel F). TSS200: 200 bp 
within transcription start site; TSS1500: 1500 bp within transcription start site; UTR: untranslated region; IstExon: the first exon of gene. 
doi:1 0.1 371 /journal.pone.00961 78.g002 



enriched (GSEA reports the gene sets with FDR <0.25) in NFAs, 
and no such gene sets were detected in FAs (Table S3, File S6). 
The most significant gene set was related to ion-channel activity 
signaling. Twenty core enrichment genes are shown in Table S4 
and the pathway heat map is shown in File S7. The genes enriched 
in the pathway encoded K + , GT~ and Ca 2+ ion-channels, which 



are known to play important roles in hormone secretion [27,28]. 
In order to determine whether there were any pivotal genes 
driving these pathways, we performed a leading edge analysis of 
the four significant gene sets. Two core genes, KCMMB4 
(potassium large conductance calcium-activated channel, subfam- 
ily M, beta member 4) and CACJVA1C (calcium channel, voltage- 
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the island, respectively), CpG island shelves (0-2 kb from shore edge, N and S indicate the upstream and downstream of the island, respectively) and 
non-CpG island probes. C and D, Relative distribution and count of the promoter associated CpGs (total 360), which were singled out from the global 
significant CpG sites (total 3027). E, KCNAB2 (cg1 81 92083) MethyLight assay in NFA and FA, which was consistent with the genome-wide DNA 
methylation analysis. 
doi:1 0.1 371 /journal.pone.00961 78.g003 
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dependent, L type, alpha 1C subunit) overlapped in the significant 
gene sets, and the ion-channel activity signal gene set highly 
overlapped with the substrate-specific channel activity signal gene 
set (Figure A and B in File S8). These data suggest that ion- 
channel pathways may play an important role in the process of 
hormonal secretion in functional PAs. As such, identified 
hypermethylation pathway genes, and in particular KCNAB2, 
may provide new insight to the epigenetic differentiation of FAs 
and NFAs. 

Association of DNA methylation with gene expression in 
sporadic PAs 

Gene expression can be modulated in multiple ways, from 
transcription to post-translational modification. Epigenetic mod- 
ification may affect gene expression levels at multiple steps, 
including DNA methylation, histone modification, chromatin 
remodeling and transcription factor interaction. To examine 
whether DNA methylation was associated with gene expression in 
PAs, we performed differential gene expression analysis among the 
24 subjects using RNA-Seq. 

Comparison between NFAs and FAs showed that gene 
expression was not found to be significandy different (FDR< 
0.01) between the two groups. However, when relaxing the 
threshold to claim statistical significance, 27 genes (FDR<0.05) 
displayed a potentially different expression profile (Table 3). Based 
on the globally significant methylated CpG sites, we selected the 
relative hypermethylation genes in NFAs and investigated their 
gene expression patterns. We found that the trend of gene 
expression was negatively associated with DNA methylation 
(Figure 4A). When matching those 27 genes with significant genes 
identified by DNA methylation analysis, three genes, including 
ribosomal protein S6 kinase, 90 kDa, polypeptide 2 (RPS6KA2), 
retinol dehydrogenase 10 (RKH10), and odontogenic ameloblast 
associated protein (ODAM) were identified with altered methyla- 
tion CpG sites in the gene body (Figure 4B). We found that genes 
RDH10 and ODAM, which were hypomethylated in the CpG 
island and enhancer region, showed increased expression in FAs 
compared to NFAs. However, the RPS6KA2 gene, which was 
hypermethylated in a CpG site in the gene body, also showed 
increased expression in FAs compared to NFAs. Finally, we found 
that ODAM was expressed in 16.7% (1/6) of FAs but showed no 
expression in NFAs (0/18) (Figure 4C). These results indicate that 
the function of DNA methylation seems to vary with context, and 
the relationship between DNA methylation and transcription is 
more complex than previously expected [29] . 

Discussion 

Genome-scale DNA methylation screening of PAs was per- 
formed to investigate whether DNA methylation was associated 
with PA invasion and histopathological subtype classification. To 
our knowledge, this is the first genome-scale DNA methylation 
analysis of sporadic PAs, demonstrating that epigenetic modifica- 
tion of key gene substrates may in part account for functional 
differentiation of PAs, and that DNA methylation analysis of key 
candidate genes may potentially be used to complement PA 
current histopathology subtype classification systems. 

Tumor invasion has been observed to occur in up to 85% of 
surgically-resected PAs, based on microscopic analysis of dural 
samples [30] . Consequently, tumor invasion is perhaps the greatest 
barrier to achieving adequate tumor control in PAs, as complete 
surgical resection of noninvasive PAs is typically achieved in 70- 
95% of noninvasive PA cases, compared with only 20^-0 % of 
invasive PAs [3 1] . In our study, no significant global differences in 



CpG methylation were identified between invasive and noninva- 
sive PAs. Although CpG island hypermethylation of a series of 
well-characterized cell cycle regulation genes, including retino- 
blastoma 1 (RBI) [32], CDKN2A [33,34], and GADD45G [35] have 
been linked to gene repression and negative regulated cell growth 
in PA, these genes were not demonstrated to contribute to the 
invasive PA phenotype in our study. Importandy, no significant 
differences in CDKN2A methylation were observed between 
noninvasive and invasive NFAs in our study, a finding which is 
consistent with prior observations [7], and may also suggest that 
CDKN2A inactivation is alternatively related to PA subtype and 
size [36]. The current findings suggest that genome-scale DNA 
methylation assessment may be utilized to identify candidate genes 
involved in PA invasion that would otherwise possibly be 
concealed in targeted gene studies. In our study, a secondary 
analysis showed that the angiogenic gene SLIT3 [37] and 
oncogene FLT1 [38] may be two candidate biomarkers for 
invasive PAs, if validated in future studies with a larger PA sample 
size. Meanwhile, although a previous study showed that promoter 
hypermethylation of CDH13 and CDH1 was detected in PAs but 
not in normal pituitary tissue, and promoter hypermethylation of 
CDH13 was observed more frequendy in invasive PAs [1 1], there 
was no significant global methylation difference between invasive 
and noninvasive PAs in our study. We did, however, find gene 
body hypermethylation in CDH7, a type 2 classical cadherin from 
the cadherin superfamily, when comparing noninvasive to invasive 
NFAs (FDR = 0.01). Further validation in a larger PA sample 
population is nevertheless warranted to identify and confirm 
genomic and epigenomic pathways involved in tumor invasion. 
Meanwhile, other epigenetic mechanisms (i.e. histone acetylation) 
should be explored in an effort to explain the character of PA 
invasion and prioritize genes associated with this clinically 
significant phenotype. 

The anterior pituitary gland (adenohypophysis) secretes six 
known hormones that regulate homeostasis, including adrenocor- 
ticotropic hormone (ACTH), growth hormone (GH), prolactin, 
thyroid-stimulating hormone (TSH), follicle-stimulating hormone 
(FSH) and luteinizing hormone (LH) [39]. Clinically nonfunction- 
ing PAs are not able to synthesize and/or secrete functional 
hormones, although in a majority of cases they nevertheless 
demonstrate positive immunostaining for LH, FSH, and/or the 
alpha-subunit [40]. PAs are monoclonal in origin and typically 
benign tumors, suggesting they arise from expansion of single 
precursor cells that possess a unique proliferative advantage [41]. 
These monoclonal adenomas therefore secrete specific hormones 
reflective of their differentiated cell of origin. In the current study, 
hierarchical clustering analysis readily separated somatotroph 
(GH) adenomas from nonfunctional adenomas, suggesting a 
potential role of DNA methylation in the differentiation and/or 
functional regulation of these tumors. Furthermore, an unexpected 
and interesting finding in the current study was that the functional 
corticotroph adenoma causing Cushing's disease and the silent 
corticotroph adenoma clustered to the same hierarchical group 
(group 2). 

The exact mechanisms by which nonfunctioning PAs retain 
immunostaining but fail to produce clinically-active hormones 
remains to be determined. Our genome-scale screening of 
variability in DNA methylation between functional and nonfunc- 
tional PAs suggests that nonfunctional PAs are globally hyper- 
methylated compared to functional ones. This finding is partic- 
ularly interesting in that the most differentially methylated gene, 
KCNAB2, encodes a potassium ion-channel that has been 
previously implicated in endocrine function pertaining to insulin 
secretion [42]. Pituitary cells resemble neurons and muscle fibers 
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Figure 4. Gene expression. A, Linear regression analysis showed negative trend of DNA methylation and gene expression. B, Three overlapped 
significant genes of DNA methylation and gene expression when compared NFAs to FAs. C, gene ODAM expression in one FA subject of d0079 which 
secreted hormones of GH and TSH, GAPDH gene expression was used as control. 
doi:10.1371/journal.pone.0096178.g004 



in that they also fire action potentials (APs) [43], which are 
mediated via expression of numerous voltage-gated sodium, 
potassium, calcium and chloride channels. These APs are 
accompanied by a rise in intracellular calcium and spontaneous 
electrical activity that drives intracellular calcium concentrations 
above the threshold for stimulus-secretion and stimulus-transcrip- 
tion coupling. KCMAB2 is a subunit of the shaker-related voltage 
dependent potassium channel, and upon binding to K + channel 
alpha subunits, contributes to regulation of channel excitability. In 
addition, KCJVAB2 has been demonstrated as an functional aldo- 
ketoreductase (AKR) [44,45], Among the primarily hypermethy- 
lated genes in NFAs, we also identified four significant genes 
PSEN2 [46], MD.N1 [47], PCSK6 and NME9 [48], which were 
relatively hypormethylated in NFAs compared to FAs (Table 2). 
Importantly, PCSK6 has been reported to be a member of the 
mammalian subtilisin-like proprotein convertase family, which 
participates in maturation of precursor protein, and is expressed at 
high levels in the anterior pituitary gland [49] . This link potentially 
explains the significant DNA methylation differences in PCSK6 
methylation between NFA and FA, and suggests that PCSK6 may 
therefore be an important biomarker for distinguishing FA from 



NFA. Collectively, these data suggest that genome-scale DNA 
methylation analysis may provide a practical, complementary 
molecular correlate to standard PA classification schemes. Further 
analysis in a larger sample size is required to support the validity of 
these findings. 

Additionally, our study demonstrates that the ion-channel 
activity signal pathway shows a relative hypermethylated profile 
in nonfunctional PA subjects. In particular, the pivotal gene 
KCMMB4 encodes an important subunit of potassium channels and 
is expressed throughout the brain (especially in the thalamus and 
the brainstem) [50]. In addition, transcription of the calcium 
channel gene C4CjV/17£" has been reported to undergo epigenetic 
regulation via DNA methylation [51], Taken together, these 
findings suggest that modulation of ion activity pathways may 
contribute to defective hormone secretion in nonfunctional PAs, 
and possibly serve as novel candidate biomarkers in PAs. 

Finally, in order to investigate the association of DNA 
methylation and gene expression, total RNA sequencing and 
differential expression analysis was performed. Although the 
expression of KCNAB2 and other genes in ion-channel activity 
signal pathway was not significantly different between FAs and 
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Table 3. List of genes with most significant differences in expression between NFAs and FAs. 



Gene 


Locus 


FA-FPKM* 


NFA-FPKM 


p value 


FDR 


ASAP2 


chr2:9346893-9545812 


41.37 


6.20 


5.0x10 ~ 5 


0.03 


ATP1B2 


chrl 7:7554253-7561 089 


9.37 


99.48 


5.0x10 ~ 5 


0.03 


CDK18 


chrl :205473683-205501 921 


3.95 


61.97 


5.0x10 ~ 5 


0.03 


CNPY1 


chr7:l 55293952-1 55326539 


0 


0.44 


5.0x10 ~ 5 


0.03 


COL9A2 


Chrl :407661 62-40782939 


91.12 


5.10 


5.0x10 ~ 5 


0.03 


CSH1 


chrl 7:61 972267-61 974021 


0.57 


0 


5.0x10 ~ 5 


0.03 


CSHL1 


chrl 7:61 986964-61 98861 8 


1.09 


0 


5.0x10 ~ 5 


0.03 


FAM138E 


chrl 5:1 02495087-1 02496558 


0 


0.49 


5.0x10 ~ 5 


0.03 


FDCSP 


chr4:71091787-71 100968 


2.05 


0 


5.0x10 ~ 5 


0.03 


CNB3 


chrl 2:6949374-6956557 


25.85 


287.08 


5.0x10 ~ 5 


0.03 


ITPR2 


chrl 2:26488284-269861 31 


2.88 


18.21 


5.0x10 ~ 5 


0.03 


MIR7-2 


chrl 5:891 55055-891 551 65 


0 


10.86 


5.0x10 ~ 5 


0.03 


ODAM 


chr4:71 062243-71 070293 


0.77 


0 


5.0x10 ~ 5 


0.03 


ONECUT1 


chrl 5:53049352-53082209 


0 


2.78 


5.0x10 ~ 5 


0.03 


PADI3 


chr1:17575592-17610727 


0.77 


0 


5.0x10 ~ 5 


0.03 


PITX2 


chr4:111538579-111563279 


16.01 


146.53 


5.0x10 ~ 5 


0.03 


PTH 


chrl 1:13513600-13517567 


0 


1.71 


5.0x10 ~ 5 


0.03 


RDH10 


chr8:74206836-74237520 


41.71 


4.72 


5.0x\10 ~ 5 


0.03 


RPS6KA2 


chr6:166822853-1 67275771 


37.64 


6.46 


5.0x10 ~ 5 


0.03 


SEC14L3 


chr22:3085521 5-30868034 


1.29 


0 


5.0x10 ~ 5 


0.03 


SEMA7A 


chrl 5:74701 629-74726299 


1.34 


83.58 


5.0x10 ~ 5 


0.03 


SLC10A4 


chr4:48485359-48491 541 


6.17 


109.53 


5.0 x 10 ~ 5 


0.03 


SMARCD3 


chr7:1 50936058-1 5097423 1 


37.70 


1 69.90 


5.0x10 ~ 5 


0.03 


SPRR2F 


chrl :1 5308461 2-1 53085989 


0 


0.60 


5.0x10 ~ 5 


0.03 


SSC5D 


chrl 9:55999869-56030466 


6.20 


81.08 


5.0x10 ~ 5 


0.03 


TFF1 


chr21 :43782390-43786644 


2.09 


0 


5.0x10 ~ 5 


0.03 


UNCX 


chr7:1272653-1276613 


0 


2.56 


5.0x10 ~ 5 


0.03 



a FPKM: fragments per kilo base of exon per million fragments mapped. 
doi:1 0.1 371 /joumal.pone.00961 78.t003 



NFAs, our global correlation analysis showed that DNA 
hypermethylation was negatively associated with gene expression. 
The elevated expression of ODAM in FAs provides a potentially 
novel insight to the pathogenesis of epithelial neoplasms, since 
ODAM is a developmental antigen with an essential role in tooth 
maturation and the pathogenesis of various epithelial neoplasms 
[52]. Although we explored differences in RNA expression of 
KCJVAB2 between NFAs and FAs, the result was not significant. A 
larger sample size would likely be needed to validate this 
association. Because gene expression is modulated not only by 
epigenome modification, low levels of gene expression cannot be 
attributed solely to hypermethylation. Additionally, the intra- 
tumoral multifocality and heterogeneity of DNA methylation of 
the PAs may account for this abberant association of DNA 
methylation and gene expression [53]. There are several 
limitations to the current study, including a relatively small sample 
size and nonuniform representation of functional to nonfunctional 
PAs and various histopathological PA subtypes. Furthermore, 
MethyLight validation of only selected candidate genes was 
performed, mandating further validation in a larger sample size 
and using another technique. The mechanisms by which 
methylation alterations of selected candidate genes result in 
observed invasion or hormonal phenotypes were also not studied, 
and remain beyond the scope of the current study. Finally, the 



current study is limited in that no true control samples (normal 
pituitary gland) were included in the analysis. Nevertheless, the 
current study suggests that DNA methylation analysis can be used 
to provide valuable insight into the phenotype of histopathological 
subtype in PAs, and if validated may complement current 
pathological classification systems. The current study is the first 
genome-wide DNA methylation analysis in PAs, and can be used 
to appropriately design and power future studies with a larger 
sample size in order to validate many of the preliminary findings 
from our study. 

In conclusion, genome-scale DNA methylation profiling and 
RNA sequencing of PAs identified DNA methylation variations in 
candidate genes associated with functional subtype (both globally 
and in selected genes) and possibly invasion. Hierarchical 
clustering analysis showed PA clustering according to functional 
status and immunohistochemical subtype, suggesting that DNA 
methylation analysis may possibly provide a clinically useful and 
complementary molecular correlate to standard PA classification. 
Differential methylation of cell motility related genes (such as 
FLT1 and SLITS) require further validation prior to being 
considered candidate biomarkers for PA invasion. DNA hyper- 
methylation of KCJVAB2 and enrichment of DNA methylation in 
ion-channel activity signal pathways may be associated with the 
endocrine-inactive status of nonfunctional PAs. 
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Supporting Information 

File SI Analysis strategies. Figure A, Filtering strategy for 
Illumina HM450 Probes. SNP10: SNPs within 10 bp from CpG 
site; SNP15: SNPs within 15 bp from the CpG site lying entirely 
within a repeat region; GhrX: chromosome X. Figure B, 
Additional filters were used in the secondary analysis of the 
invasive phenotype of NFAs. A total of 884 probes were enrolled 
in further analysis, of which 355 were enhancer (Enh)- associated, 
20 were promoter (Pro)- associated, and 4 were both (Bot) 
promoter- and enhancer-associated, and we identified 34 signif- 
icant (Sig) CpGs that were independently associated with 
enhancers and were hypomethylated in invasive NFAs compared 
to noninvasive NFAs. 
(TIF) 

File S2 The distribution of probes on the Illumina 
HM450 DNA methylation platform. Figure A, The distribu- 
tion of 383,718 HM450 probes stratified by CpG density (in and 
out of CpG islands). Figure B, The distribution of differentially 
methylated CpGs across different gene-related and intergenic 
regions. 
(TIF) 

File S3 Hierarchical Clustering of DNA methylation in 

24 PA cases. The HM450 probes with the highest 1% of 
standard deviation across the set of 383,7 18 global HM450 probes 
were used. 
(TIF) 

File S4 Hierarchical Clustering of DNA methylation in 

24 PA cases. The HM450 probes with the highest 5% of 
standard deviation across the set of 383,7 18 global HM450 probes 
were used. 
(TIF) 

File S5 Hierarchical clustering of DNA methylation 
levels within the KCNAB2 gene locus. Beta values with 
50% of the highest standard deviation were used. Pro-a: promoter 
associated. Compared with FAs, NFAs showed relative hyper- 
methylation across nearly the whole gene, especially in promoter- 
associated CpG sets. 
(TIF) 

File S6 Enrichment plots of significant gene sets. The 

score at the peak of the plot is the enrichment score for the gene 

set. "0" represents NFAs, and "1" represents FAs. 

(TIF) 

File S7 Heat map of genes in the ion-channel activity 
signal gene set. DNA methylation values are represented as 
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colors, with red representing DNA hypermethylation and blue 

representing DNA hypomethylation. 

(TIF) 

File S8 Leading edge analysis of the four overlapping 
significant gene sets. Figure A, Darker color represents greater 
overlap between the subsets. The gene set ion-channel activity 
signal pathway had a lot of overlap with the substrate specific 
channel activity signal. Figure B, KCNMB4 and CACMA1C were 
enriched in all of the four gene sets. DNA methylation values were 
represented as colors, where the range of colors (red, pink, light 
blue and dark blue) shown the range of methylation values (high, 
moderate, low and lowest). 
(TIF) 

Table SI Top 17 significant genes with most differen- 
tially-methylated enhancers between invasive and non- 
invasive NFAs. 

(XLSX) 

Table S2 Hypermethylation genes in NFAs, which are 
consistent with results of Illumina Infinium Methylation 
27K arrays. 

(XLSX) 

Table S3 Significant gene set enrichment in NFAs. Four 
significant pathways were listed in the table, and ion-channel 
activity pathway showed high enrichment score with the lowest 
FDR value. 
(XLSX) 

Table S4 Gene subset of ion-channel activity signal 
pathway. The majority of genes were associated with potassium 
channels. 
(XLSX) 

Table S5 Full names of genes mentioned in the text. 

(XLSX) 
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